Backward error analysis and the substitution law for Lie 

group integrators 

Alexander Lundervold * Hans Munthe-Kaas t 



Abstract 

Butcher series are combinatorial devices used in the study of numerical methods for differ- 
ential equations evolving on vector spaces. More precisely, they are formal series developments 
of differential operators indexed over rooted trees, and can be used to represent a large class 
of numerical methods. The theory of backward error analysis for differential equations has 
a particularly nice description when applied to methods represented by Butcher series. For 
the study of differential equations evolving on more general manifolds, a generalization of 
Butcher series has been introduced, called Lie-Butcher series. This paper presents the theory 
of backward error analysis for methods based on Lie-Butcher series. 

1 Introduction 

A fundamental tool in the field of numerical integration of ordinary differential equations on R n 
is the theory of Butcher series (B-series). These are formal series expansions of vector fields 
and flows, expanded over the set of rooted trees. Many numerical methods can be formulated in 
terms of B-series, and they can be used to, for example, study order theory, structure preserving 
properties of integrators, backward error analysis and modified vector fields [3l [T7 l 151 I7[ [TO ] \9\ [23] . 
In the more general setting of differential equations of the form 

y' = F(y), yeM, F : M -> TM, (1) 

where M is a (homogeneous) manifold and F a vector field on M, the role of B-series is played by 
the Lie-Butcher series (LB-series) [2E1 1221 123] • Considering the importance of classical B-series, 
LB-series are objects of great interest. 

The B-series are based on the elementary differentials associated to vector fields, and these 
can be constructed as homomorphisms from the free pre-Lie algebra (or Vinberg algebra) into the 
pre-Lie algebra of vector fields [5J. In the setting of LB-series we get a similar picture, only now 
the pre-Lie algebras are replaced by the so-called post-Lie algebras, defined in [Ml 12"?] . 

In the present paper we will explore the substitution law for Lie-Butcher series, formulated 
in the language of enveloping algebras of post-Lie algebras: the D-algebras of [35]. Once the 
substitution law is understood, it can be applied to backward error analysis. The basic idea of 
backward error analysis is to interpret the numerical solution of a differential equation as the 
exact solution of a modified equation, and then use this equation to study the numerical method. 
Analogous to classical backward error analysis (as developed in [THl 03 [5J \5\ ) , its generalization to 
the Lie group setting has a particularly nice description for methods based on Lie-Butcher series. 

Note that the construction of series expansions in the present paper is purely formal: there will 
be no study of convergence. This separation between the algebraic and the analytic framework for 
backward error analysis is also present in the setting of B-series, where the main algebraic references 
are [T71 [5J [51 [5] and the analytic references are [H [321 H3 • An analytic study of backward error 
analysis for Lie group methods can be found in [13] , 
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The present study of the backward error and substitution law for Lie group integrators is 
interesting from a purely algebraic point of view, as this work provides an explicit description 
of automorphisms of post-Lie algebras. From a numerical point of view, the theory has several 
applications. The algebraic structures of backward error analysis is important in the analysis 
of numerical integration algorithms. Additionally, in the case of classical B-series, such algebraic 
techniques have recently been applied more directly as a computational tool 8] . Similar techniques 
in the setting of Lie group integrators is a promising approach to structure preserving integration 
of problems of computational mechanics, such as Lie-Poisson systems. 

2 Lie— Butcher series 

In this section we will define D-algebras, and show how they give rise to Lie-Butcher series. In the 
next section we will apply them to the study of the substitution law and backward error analysis 
for Lie group integrators on manifolds. 

2.1 Trees and D-algebras 

Ordered rooted trees and forests. Some basic definitions follow. For a more comprehensive 
introduction to the combinatorics of trees applied to numerical integration, see [4] or [19]. Let OT 
denote the alphabet of all ordered (i.e. planar) rooted trees: 



The root is the bottom vertex and we consider the trees to grow upwards from the root. The 

trees being ordered implies that . This is different from classical B-series theory, where 

the order of the branches is of no significance. Let OF denote the set of ordered forests, i.e. all 
possible empty and non-empty words written with letters from the alphabet OT: 



where I denotes the empty word. On OF we define the concatenation product u>%,u>2 ^ w i w 2 j which 
creates a longer word by joining uj\ and UJ2 end-to-end. This is an associative, non-commutative 
product with unit I. Let B + : OF — > OT denote the operation of adding a root to a word, e.g. 

All of OF is generated from I by concatenation and adding roots. The order 
of a forest, \u>\ — \ri...T k \, is defined by the recursion |I| = 0, \ri...r k \ = \t\\ + ■•• + \r k \, 
\B + oj\ = \uj\ + 1, i.e. the order counts the number of vertices in a forest. Let k be a field of 
characteristic 0, e.g. k = R or k = C. The k-vector space of all finite k-linear combinations of 
elements in OF is the non-commutative polynomial ring over OTQ denoted by Af = k(OT). The 
k-vector space of infinite linear combinations of OF is Af* — k((OT)). Af* is the dual space of Af, 
with the dual pairing (■,■): Af* x Af k defined such that the words in OF form a orthonormal 
basis: (u^u^) = if u>\ ^ uj 2 , and ((jJ,lj) = 1. Thus for a £ Af* we have a(uj) = (a,uj) and 
a = X^gof a ( CJ ) w - I n the latter sum we understand Af* as the projective limit Af* = hrnA4, 
where Af k — spanjw £ OF: \uj\ < k}. An infinite a £ Af* is uniquely defined by its finite 
projections a k £ Af k for all k £ Z, where a k = X)|w|<fc a ( CJ ) w is the orthogonal projection of a onto 
the subspace Afk C Af. 

M with concatenation product can equivalently be defined as the linear space spanned by trees, V = fc{OT}, 
equipped with a tensor product. Hence M can be defined as the tensor algebra on V. However, because we need 
other tensor products later we prefer the definition via concatenation of words. 



OT = {.. 




...}. 
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Remark 2.1. In many applications it is necessary to generalize to spaces built from trees with 
colored vertices. The theory extends from the above presentation with only minor modifications. 
Let C be a (finite or infinite) set of colors. A coloring of a tree or a forest is a map from its 
vertices to C. Let OTc and OFc denote colored trees and forests. For each c € C we have the 
operation _B+ : OFc — > OTc creating a tree by adding a root of color c to a word. We identify 
C C OTc C OFc as the subset of single vertex trees. In the colored context we permit more 
general gradings | • | on OTc- We allow the assignment of arbitrary positive integer weights 

|c| e N to the single vertex trees C C OT c , extended to OF c by \t x . . . r k \ = \t\\ H + \tu\ and 

= \oj\ + \c\. The definitions of finite and infinite linear combinations of forests Afc = fc(OTc) 
and = k((OTc)) arc similar to the uni-color case. 

Definition 2.2. The left grafting product ■ rx ■ : Af <Ei Af — > Af is defined recursively as follows: 
let r G OT and w,wi,w 2 € OF. Then 

I rx to = oj 

t rx I = 

lo rx • = B + (u>), 

t rx (LO1UJ2) = (t rx uji)lu2 + u>i(t rx UJ2) 

(tu) rx lo\ = t rx (oj rx loi) — (r rx lo) rx uj\ 

The product is extended to all of M and Af* by linearity and projective limits. 
For example, 

~rxl = *V + 2*} + Y 

Note that grafting satisfies a Leibniz rule with respect to the concatenation product. If we define 
touj = tlj + t rx oj, we see that r rx (oj rx oj\) = (row) rx oj\. More generally, oj\ rx (0J2 rx oj) = 
(0J1O0J2) rx oj, where o is the associative product defined as follows: 

Definition 2.3. The Grossman-Larson product o : Af ® Af — > Af of wi,W2 € OF is defined in 
terms of the grafting product as: 

B + (0J1O0J2) — oj\ rx B + (oj2), 

and is extended by linearity. 

It is clear that if we write wilu^] for oj\ rx 0J2, we have the following structure on Af: 

Definition 2.4 ( 28J). Let A be a unital associative algebra with product f,g <—> fg, unit I and 
equipped with a non-associative composition (.)[.] : A <E> A —> A such that I[g] = g for all g € A. 
Write T>(A) for the set of all / G A such that /[■] is a derivation: 

V{A) = {/ e A I f[gh] - (f[g])h + g(f[h]) for all g, h e A}. 

Then A is called a D-algebra if for any derivation / € ~D{A) and any g *E A we have 

(i) g[f]eV(A) 

(ii) f[g[h}} = (fg)lh] + (f[g})[h}. 

The free D-algebra. We note that a morphism T : A — > A' of D-algebras is an algebra 
morphism satisfying F(D{A)) C D(A') and J"(a[6]) = J"(a)[J"(6)] for all a,b £ A. The D-algebra 
Af plays a special role: it is a universal object. 
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Proposition 2.5 (|28j). Let OT be planar trees decorated with colors C . The vector space M = 
R(OT) is a free D-algebra over C. That is, for any D-algebra A and any map v.C^t D(A) there 
exists a unique D-algebra homomorphism T v : N — >• A such that T v (c) = v(c) for all c G C. 

C c ► TV 

3! T v 

D(A) c — ► A 

We will see that based on this result we can construct elementary differentials and Lie-Butcher 
series for Lie group integrators, and also define the substitution law. To achieve this we utilize 
D-algebra structure of differential operators on manifolds . 

The D-algebra of differential operators. There is a D-algebra based on the space of vector 
fieldtQ on the manifold M. Consider the space C°°(M, g) —: g M , where g C XM is a Lie sub- 
algebra of the set of all vector fields on M. For ^> G g M and V € g, the Lie derivative V[i>] G g M 
of ^ along V defined by 

VmP)--=j t \t=o*(eMtV)(p)). (2) 

V[-] is a first order differential operator on g M , satisfying V'lTivf'] = V^/i]^ + where h G 

C°°(M, H.) is a scalar function^] The Lie derivative gives rise to differential operators of higher 
degrees through concatenation: the concatenation of V, W € g is a second-degree differential 
operator defined by VW[^] := V[W[^]]. The C°° (M, M)-module of all differential operators, 
including the ones of higher degree, and the degree zero operator spanned by the identity operator 
I, is called the universal enveloping algebra U(g) of g. We extend the structure to the space 
e°°(M, U(g)) =: U(g) M as follows: for f,g£ U(g) M f[g] G U{g) M is defined by 

f[g}(p):=(f(p)[g})(p) (3) 

and fg G U(g) M is defined as 

fg(p) ■= f(p)g(p)- (4) 

The latter operation is called the frozen composition of / and g. For two vector fields / and g 
written in terms of the standard coordinate frame {d/dxi}, the operations take the following form: 

^ = E/*»^4f (6) 

The operations ^ and Q endows the space U(g) M with the structure of a D-algebra, where 
the derivations are the vector fields in g M : 

Lemma 2.6. Let f G g M and g : h<E U(g) M . Then 



f[gh] = f[g]h + g(f[h}) 
f[g[h}} = {fg)[h]+f[g\[h\. 



Hence U(g) is a D-algebra. 



The composition of / and g as differential operators, defined by (fog)[h] := is called 

non-frozen composition^ 

2 The vector fields are interpreted as differential operators acting on functions. 
3 This is true when g M is replaced by S M for any vector space H. 

4 We note that the two operations /, g i— > f[g] and /, g i— > fog gives U(q) m the structure of a unital dipterous 
algebra (as defined in 21 ). 
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Post-Lie algebras. The theory of Lie-Butcher series can be reformulated in terms of post-Lie 
algebras. These were first studied in the setting of operads by Vallette |34j . and also by the authors 
in |27] . Our main motivation for the construction of post-Lie algebras was their relation to the 
D-algebras defined above, which are universal enveloping algebras of post-Lie algebras. 



2.2 Lie— Butcher series 

Classical B-series. Recall (see e.g. [T7|) that a B-series is a (formal) series indexed over the 
set NT of non-planar rooted trees (i.e. trees without any ordering of the branches) and can for a 
vector field / be written as 

B hJ (a)(y) = a(I)y + £ ^\a(r)F f (r)(y). (7) 

t£NT v ; 

Here ct(t) is the symmetry factor for r £ NT, and a is a map a : NT — > R. The map Ff(r) : K™ — > 
K™ is the elementary differential of the tree r, obtained recursively from / and its derivatives: 

F f (.)(y) = f(y), F f (r)(y) = /^G/)(^/(n)(j/), ...,7>(r ro )(y)), (8) 

where r = £? + (ti, . . . r TO ) and /( m ) is the rath derivative of the vector field. The parameter h 
represents the step-size of the numerical method giving rise to the B-series. 



LB-series. We will consider a more general setting: that of differential equations evolving on 
manifolds. Let M be a manifold and XM the Lie algebra of vector fields F : M — > TM on M. 
The fundamental assumption for numerical Lie group integrators is the existence of a frame on 
TM, defined as a finite number of vector fields {-Ei, £2, ■ ■ • , -Em} spanning the tangent space T p M 
at each point p £ M. The frame is allowed to be overdetermined. It generates a Lie algebra g, 
and it is assumed that flows of vector fields in g can be computed exactly [251 ISj • Any vector 
field F : M — > TM can be written as = J2 a i^i(p)- We will study vector fields of the form 
F(p) — S fi{p)Ei(p) where fa : Af — > K arc smooth functions. Given such a vector field, let 
/ <= M be defined as f p = Y^, fi(p)Ei. We say that f p has coefficients frozen relative to the frame. 
In other words, to each such F £ XM there is an associated / £ g M so that F(p) = f p ■ p, where 
fp ■ p denotes evaluation of f p in p. We will often refer to such / £ g M as vector fields. 
The general differential equation (nj can now be written as 



The Lie-Butcher series are expansions over OT associated to integrators of this equation, just as 
B-series are associated to differential equations expressing the flow of vector fields in R™. The 
non-commutativity of combining vector fields is reflected in the planarity of the trees in OT. 

Now we can construct the elementary differentials needed to define Lie-Butcher series. As in 
the classical case they can be expressed recursively by a function T based on trees. 



Definition 2.7. The elementary differentials associated to a vect or field / 
algebra morphism Ff : Af — > U(g) AI we get from Proposition 
(i.e. C = {•} and v : • H> / in Proposition 2.5). Hence Ff is defined by 



M -> is the D- 
by associating the tree • to / 



(i) 7>(I)=I 

(ii) F f {B+{u J ))=Ff{u J )[f] 

(iii) F f (uj 1 uj2)=Ff(uj 1 )Ff(uj2) 



y' = fy-V, where / £ g 



(9) 



When the vector field / is clear from the context we will occasionally write F instead of Ff . 
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Definition 2.8. For an infinite series a € TV* = R((OT)) a Lie-Butcher series is a formal series 
in U(g) M defined as 

B f (a)= J2 ol(w)T(u). 
weOF 

For a vector field / this can also be written as the commutative diagram 

{•} c ► Af* 



M 



u(b) 



Bf 
M 



where Bf is the unique D-algebra homomorphism given by Proposition 2.5 



Remark 2.9. By coloring the vertices of the trees vi a map v we can define T and B for multiple 
vector fields. The elementary differentials T v are still obtained from Proposition |2.5[ but the set 
C will contain multiple colors. 



2.3 Some algebraic constructions 

Before we show how LB-series can be used to represent flows of vector fields on manifolds we must 
conduct a closer study of the space where the coefficients a live. To understand the various ways 
we can represent such flows it will also be helpful to look at some Lie idempotents, namely the 
eulerian and Dynkin idempotents (Section 2.3.2), and also certain non-commutative polynomials 
called Bell polynomials (Section |2.3.3 ) . We will follow the presentation in [55] and [55] . 



2.3.1 The Hopf algebras U Sh and H N 

It is well known that inserting a B-series Bhj(a) into another series Bhj(b) results in a B-series 
Bhj(a)(Bhj(b)(y)) = Bhj(a ■ b)(y). The product a ■ b on the set of maps a : OT — > K with 
a(I) = 1 gives rise to a group, called the Butcher group O [18] . This is the group of characters in a 
variant of the Connes-Kreimer Hopf algebra of renormalization [TT] [2] . A similar result holds for 
LB-series, where the Hopf algebra of Connes-Kreimer is replaced by a more general Hopf algebraic 
structure on the set of rooted trees. This Hopf algebra was introduced in [55J. See also [521 153"] . 

Note first that the vector space R(OT) spanned by trees can be turned into a Hopf algebra by 
using concatenation as product and deshuffling as coproduct. The deshuffling coproduct Ash is 
results from requiring the trees to be primitive, and extending by concatenation: 

Ash(r) = r®I + I®r, ^sh(riT 2 ) = A s/i (ti)A 5/i (t 2 ), 

where t, ti, t 2 arc trees. The antipodc is defined by 

S^TiTa • • • T„) = (-l) n T n T n -l ■ • ■ Tl, 

and the unit r] and counit e by 77(1) = I, and e(I) = 1, e(cj) = 0, for all uj e OF\I. 

The vector space Af = K(OT) can also be turned into an algebra using the shuffle product 
LU : Af <£) Af — >• Af defined recursively byIl_Uu; = k; = cjLLlI and 

(tiLJi) LU (r 2 W 2 ) = T 1 (UJ 1 LU T 2 Ld 2 ) + t 2 (t 1 uj 1 LU u 2 ) (10) 

for ti,t 2 G OT, u)\,u) 2 £ OF{^] This algebra can be given the structure of a bialgebra in several 
different ways. We can equip it with the coproduct A c : Af — > Af <S> Af given by deconcatenation of 
words: 

n-i 

A c (w) = I ® w + w ® I + ^Vi • • -n ® T l+ i • • • T n , (11) 

i=i 

5 Coproducts will occasionally be written using the Sweedler notation A(a>) = ® w (2)- 



G 



where to = t\ ■ ■ ■ r„. This results in the shuffle bialgebra, which equipped with the same antipode, 
unit and counit as the deshuffle Hopf algebra defines the shuffle Hopf algebra Hsh^\ 

If we instead equip N with the coproduct An : TV — > N^N defined recursively as Am (I) = I<8>I 
and 

Ajv(W) =u>t®I+ Ajv(o;)lu • (/ ® B+)Ajv(.B _ (t)), (12) 

where r g OT, cj € OF, we get another bialgebra Hat. Here LU • : N m -> TV (g) TV denotes shuffle 
on the left and concatenation on the right: (wi <£> w 2 ) LU • (w 3 <E> u> 4 ) = (uji LU cj 3 ) (g) (u^^)- An 
explicit description of the coproduct in terms of tree cuts can be found in Section [XT] below, and 
in [21], where it was shown that A at is the dual of the Grossman-Larson product and that Hn 
forms a Hopf algebraj^] This is the Hopf algebra governing composition of LB-series (Theorem 
2~T0l ). 

To simplify the expressions we introduce a magmatic structure (i.e. the structure of a set 
equipped with a closed binary operation with no further relations) on OF. Additional details and 
motivation for the introduction of this structure can be found in Section [4j Let ui , u>2 be two 
elements of OF and define the operation x : OF x OF — > OF by 

uji x ll>2 = ujiB + (0J2) ■ (13) 

For example, 

Vx.i=W 

This operation is magmatic, and the empty tree I freely generates all of OF. The operation is 
extended to M = R(OT) via linearity. If u) = v\ x 1*2, w^I, then we call v\ the left part ujl of to 
and V2 the right part lor. The shuffle of two elements of the magma can be defined as 

V LU UJ = (v LU LOl) X Vr + (vl LU w) X UJ R , (14) 

wLUl = lLUcj = tj, and we notice that the coproduct in Hjv can be written as 

A N (cu) = uj(g)l + A n (lu l )lu xAa^wr), (15) 
where UJx now denotes shuffle on the left, magma operation on the right. 



Characters and the composition of LB-series. Recall that a character of a Hopf algebra 
(H, A, •) over a field k is an algebra morphism a : H — > k, e.g. a{a - b) — a(a)a(b), and a(l#) = 1, 
where a, b G H , and 1h denotes the unit. The convolution product a * /3 of two characters is 
defined by 

jjAj/®jj°1k®k^k. 

This gives the set of characters G(H) of H the structure of a group. In fact, the field k can 
be replaced by any commutative algebra A, giving rise to A-valued characters. Another type of 
character we will need later are the infinitesimal characters. An A-valued infinitesimal character 
is a linear map a : H — > A satisfying 

a(h ■ h!) = n A (a{h),6(ti)) + n A (S(h),a(h')), (16) 

where ha is the product in A and 8 is the composition of the counit of H and the unit of A, 
8 = t\a ° £■ The characters and infinitesimal characters are connected via the exponential and 
logarithm, see e.g. [2~4"] . 

The group structure of the characters in "Hjv exactly corresponds to the composition of LB- 
series. 

Theorem 2.10 ([28 ). The composition of two LB-series is again a LB-series: 

Bf(a)[B f (p)]=Bf(a*l3), 
where * is the convolution product inHtq . 

6 Note that the concatenation dcshufning Hopf algebra is dual to the shuffle deconcatenation Hopf algebra. 
7 Being a graded and connected bialgebra "Hat is automatically a Hopf algebra. A more direct argument, and 
formulas for the antipode, can be found in |28| 
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2.3.2 Lie idempotents 

A Lie polynomial over an algebra A is an element of the smallest submodule of M.{A) that is closed 
under the bracket [P, Q] := PQ — QP in M.(A) . The Lie algebra of these polynomials is the free Lie 
algebra Lie(A) on A |33j . There are several important idempotent maps, called Lie idempotents, 
from R(A) to Lie(A). 

Eulerian idempotent Let H be a commutative, connected and graded Hopf algebra. Consider 
Endfc(Zf) = Homk{H, H) equipped with the convolution product *. Let Id £ Endfc(H) be the 
identity cndomorphism and 5 — rj o e £ End^ (H) the unit of convolution. 

Definition 2.11 ([20 ). The Eulerian idempotent e £ End(-ff) is given by the formal power series 

t*2 t*3 T*i 

e := log* (Id) = J -_++.. . — + -.., 

2 6 l 

where J = Id — 5. 

Proposition 2.12 (|20]L For any commutative graded Hopf algebra H , the element e £ Endj.(iJ) 
defined above is a Lie idempotent in H . That is, eoe — e and it has image in the free Lie algebra. 

The practical importance of the Eulerian idempotent in numerical analysis arises in the study 
of backward error analysis, where the following lemma provides a computational formula for the 
logarithm: 

Proposition 2.13 f[2"2"]). For a £ G(H) and h £ H, we have 

log*(a)(h)=a(e(h)). 

In other words, the logarithm can be written as right composition with the eulerian idempotent: 

log* =_oe:G{H)-+ S (H). 

Dynkin idempotent The classical Dynkin operator on the shuffle Hopf algebra is given by 
left-to-right bracketing: 

D(ai...o n ) = [••• [[01,02], 03], . . . ,o„], where [ai,a 3 ] = a^a,- - Ojflj. 

Letting Y(u>) = #(cj)o; denote the grading operator, where is word length, it is known 

that the Dynkin idempotent Y^ 1 o D is an idempotent projection on Lie(A). As in |12j . the 
Dynkin operator can be written as the convolution of the antipode S and the grading operator Y: 
D = S *Y . This description can be generalized to any graded, connected and commutative Hopf 
algebra H: 

Definition 2.14. Let H be a graded, commutative and connected Hopf algebra with grading 
operator Y : H — > H. The Dynkin operator is the map D : H — > H given as 

D := S * Y. 

2.3.3 The non-commutative Bell polynomials 

In [2 6) some non-commutative polynomials B n were introduced to express the Butcher order theory 
of Runge-Kutta methods on manifolds. In [35] it was observed that these polynomials were a non- 
commutative analogue of Bell polynomials, and that they could be used to study more general 
flows on manifolds. We recall their definition here. 

Let I = {dj}°^ 1 be an infinite alphabet in 1-1 correspondence with N + , and consider the free 
associative algebra V = R(Z) with the grading given by \dj\ — j and \dj 1 ■ ■ ■ dj k \ = ji + ■ ■ ■ + jk- 
Let d: T> — > T> be the derivation given by d{di) = di + \, linearity and the Leibniz rule d{uj\UJ2) = 
d(oJi)uj2 + U)\d{u)2) for all wi, w 2 £l* . We let #(w) denote the length of the word ui. 
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Definition 2.15. The non-commutative Bell polynomials B n := B n {d\, . . . ,d n ) £ T> are defined 
by the recursion 

Bo = I 

B n = (di + d)B n _ x = (di + d) n l for n > 0. 

The first few are 

Bo = I 

Bi = di 

B 2 = d\ + d 2 

B 3 = d? + 2did 2 + d 2 di + d 3 

B 4 = d\ + 3d^d 2 + 2did 2 d 1 + d 2 d\ + 3did 3 + d z d x + 3d 2 d 2 + d 4 . 

We write -B„ & := B n j : (d 1 , . . . , d n _fc + i) for the part of B n consisting of the words of length k, e.g. 
^4,3 = 3dfd 2 + 2d\d 2 di + d 2 d\. It is often useful to employ the polynomials Q n and Q n ^ related 
to B n and by the following rescaling: 

Q„,fc(di,...,d n _fe + i) = -^.B,, ,fe(l!di, . . . , jldj, . . .) = V" k(cj)uj (17) 

M=n,#(w)=fc 

Qn{di 7 . . . , d n ) = Qn,k{d\, ■ ■ ■ , d n _fc + i) (18) 

fc=l 

Qo := I. (19) 

These polynomials can be used to define an operator Q on any graded Hopf algebra H. Let d^ be 
defined on H* by 

di(a) = on = a\ H . , didj(a) — a t * a 3 , (20) 

where a, is the degree z component of a and * is the convolution product. The operator Q is a 
bijection from infinitesimal characters to characters of H (for details, see [2"2"]). 

2.4 Lie— Butcher series and flows of vector fields 

Flows i/o ^ y(t) = ^t(yo) on the manifold M can be represented by LB-series in several different 
ways. Here are three procedures, giving rise to what we will call LB-series of Type 1, 2 and 3: 

1. In terms of pullback series: Find a <= G(Hn) such that 

*(y(t)) =B(a)(lto)[*] for any * e C/( fl ) M . (21) 

This representation is used in the analysis of Crouch-Grossman methods by Owren and 
Marthinsen [31]. In the classical setting this is called a S'-series [2"5] . 

2. In terms of an autonomous differential equation: Find /3 £ g(Hiy) such that y(t) solves 

*/(*)= S(/3)0/(i)). (22) 



This is called backward error analysis (confer Section 3.3). 

3. In terms of a non-autonomous equation of Lie type (time dependent frozen vector field): 
Find 7 6 Q(Hsh) such that y(t) solves 

y'(t) = (jLBfrXvo)) y(t). (23) 

This representation is used in [S5J • In the classical setting this is (close to) the standard 
definition of i?-series. 
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The algebraic relationships between the coefficients a, (3 and 7 in the above LB-series are [22] 



P = aoe e is Euler idempotent in Hn- (Proposition 2.131 

a = exp°(/3) Exponential wrt. GL-product 

7 = aoY~ l oD Dynkin idempotent in Hsh- [HI Proposition 4.4] 



a = Q(j) Q-operator (20) in H S h- [H Proposition 4.9] 



By using these relationships one can convert between the various representations of flows. 

In the notation in the following examples of LB-series we suppress the vector fields and ele- 
mentary differentials, and phrase the LB-series in terms of the (dual of the) coefficient functions. 

Example 2.16 (The exact solution). The exact solution of a differential equation 

V'(t) = F(y(t)) 

can be written as the solution of 

y' = Ft-y, y(0)=y o , 

where F t — F(y(t)) S Q is the pullback of F along the time dependent flow of F. Let F t — ^Bt("f). 
By [221 Proposition 4.9] the pullback is given by B t (Q{jExact))[F], so 

^°7Exact = <3(7Exact)[»] => 7Exact = Y~ 1 oB + (Q(7Exact ) ) ■ 

Note that this is reminiscent of a so-called combinatorial Dyson-Schwinger equation |15j . Solving 
by iteration yields 




^r + V + 2^ + Y + !] +i(^ + V 



7Exact — 



Note that a formula for the LB-series for the exact solution was given in [31] . We observe that 
there cannot be any commutators of trees in this expression. Therefore, in LB-series of numerical 
integrators, commutators of trees must be zero up to the order of the method. 

Example 2.17 (The exponential Euler method). The exponential Euler method [T!5] can be 
written as follows: 

y n+ i = cxp(hf(y n ))y n , 

or, by rescaling the vector field /, as 

Un+i = exp(/(y n ))y n . 

This equation can be interpreted as a pullback equation of the form $(y„ + i) = £>(exp(»))[$]y„, so 

, . „ 1 1 
a = exp(«) = I + • + + ^««« H 

(Here the Grossman-Larson product is the same as concatenation). Note that exp(«) = Q(«), so 
the Type 3 LB-series for the Euler method is simply 

7Eulcr = •• 
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Example 2.18 (The implicit midpoint method). The implicit midpoint method [19] can be pre- 
sented as: 



v = /(exp(^o-)y„) 
y n+ i = cxp(er)y„ 

We make the following ansatz: 

CT = ^aH^ = a(.). + a(T)I + a(Y) Y + a ([•,!]) [•,!}+ a !+••• 

i.e. that a can be written as an infinitesimal LB-series. From Equation [24j we get that 

v ^ 3 r 1 

3=0 J 



(24) 
(25) 



(26) 



(27) 



Since there are no forests in this expression, we must have a([ui, u)']) = for all uj, lj' G OT. If we 
write t = B + (ti • ■ ■ Tj), then by combining Equation 27 with the ansatz, we see that coefficients 
of the LB-series are given recursively as «(•) = h, 



(28) 



Hence 



^Midpoint — 



it i i 



2! 2 4 




3 Substitution law for Lie— Butcher-series 

In this section we will generalize the substitution law for B-series [7] to LB-series. Once the 
substitution law has been established we will apply it to backward error analysis for numerical 
methods based on LB-series. 



3.1 The substitution law 

Consider AT as a D-algebra where the derivations are the Lie polynomials D(j\f) = q{iHn) HJV. 
By the universal property of Af, we know that for any map a: C — > D(j\f) there exists a unique 
D-algebra homomorphism T a : M — > M such that F a (c) = a(c) for all a € C. This is called the 
substitution law. 



Definition 3.1. For any map a : C — > D(Af) the 
unique D-algebra homomorphism a* : j\f — > j\f 
such that a(c) = a * c for all c g C is called 
a-substitutior|3 



C c ► Af 



D(AT) 



AT. 



Theorem 3.2. The substitution law defined in Definition \3.1\ corresponds to the substitution of 
B-series in the sense that 

B Bf -0J)(a) =B f (p*a) 



s In most applications we want to substitute infinite series and extend a* to a homomorphism a* : jV* — > jV* . 
The extension to infinite substitution is straightforward because of the grading, we omit details. We write a* also 
for infinite substitution. 



11 



The theorem is easily proven by using the following lemma: 



Lemma 3.3. For all (3 : {•} -> D(Af*) and all B-series Bf : Af* -> U(g) M , the composition 
Bf o (3 has image in q m . In other words, B-series maps D(Af) to derivations on M. 

Proof. It is enough to prove this for Lie polynomials [^j Since Bf is a D-algebra homomorphism 
it maps trees to derivations, so the only thing we have to check is that the commutator [V, W] = 
VW — WV of two derivations V and W is a derivation. This is a straightforward calculation. □ 



Proof of Theorem 3.2 Except for the use of Lemma 3.3 the proof is purely categorical. Let Bf 
be a B-series. The composition of Bf with the map /3-k can be written in diagrammatic form as 



{•} 



Af* 



D{Af) 



/3* 



Af* 



U(Q) 



Bf 

M 



By Lemma 3.3 the composition of the two diagonal arrows and Bt actually has image in q 



Therefore the universal property for the diagram obtained by adding the map Bf o 
to the above diagram shows that Bf o /3-k = Bs f0 p, and hence the theorem. 



M 

□ 



Many of the useful properties of the substitution law follow immediately from the fact that a* 
is a D-algebra homomorphism. For example, a* : Af — > Af is a linear map which for any n, nf £ Af 
satisfies 



a -k (nn!) 
a-k (n rx n) 
a -k (non ) 
a -k (noS) 
a * (noe) 



(a -kn)(a-k n') 
(a-k n) rx (a* n) 
(a * n)o(a 7k- n ) 
(a * n)oS 
(a -k n)oe 



where S is the antipode and e is Euler map in %n- 

The free D-algebra Af is the universal enveloping algebra of the free post-Lie algebra q of rooted 
trees [27] . By defining a coproduct by requiring that the elements of q are primitive (e.g. the 
deshuffle coproduct of Section 2.3.1), it is a bialgcbra. The unique D-algebra morphism a-k is a 
coalgebra morphism for this coproduct: 

Lemma 3.4. The map a-k is a coalgebra morphism with respect to the coproduct given by deshuf- 
fling of words (Section 2.3.1). That is, 

(a -k tgia-k) o Ash — Ash ° a*, 

where Ash denotes the deshuffling coproduct. 

Proof. The result is easily proven for primitive elements. The general case follows by induction 
on the length of words. □ 



'Lie series are formal series whose homogeneous components are Lie polynomials |33| 
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Remark 3.5 (The Hopf algebra for the substitution law). Based on the results in [5] and the fact 
that the operad governing post-Lie algebras is known, it is possible to describe the Hopf algebra for 
the substitution law following the program in This is a project currently under development 

3.2 A formula for the substitution law 

The substitution law can be calculated recursively using a formula involving trees. To write down 
the formula we need to look at cutting operations on trees and forests. 

Cutting trees and forest. Let r € OT be an ordered rooted tree. An elementary left cut c of 
r is a choice of a set of branches E of r to be removed from r. These are chosen in a systematic 
manner: if an edge e is in E then all the branches on the same level and to the left of e must 
also be in E. Each cut splits r into two components: the pruned part P^(r) consisting of the 
trees that were cut off concatenated together, and the remaining part i?g Z (r) consisting of the 
tree containing the root. We also consider the empty cut, i.e. the cut c so that P^{t) = I and 
R c el {r) = r, to be an elementary cut. 



r-4 








I 


I 




• 


• 




• 


V 



A left admissible cut on r consists of a collection of elementary cuts applied to r with the property 
that any path from the root to any vertex of r crosses at most one elementary cut. The pruned 
parts corresponding to each elementary cut are shuffled together, with no internal shuffling of the 
trees resulting from each elementary cut. An admissible cut of a tree results in a collection of 
shuffles of forests P c {t) and a tree R c {t). The collection of all left admissible cuts for a tree r is 
written as LAC(t). 



T-4 




R c (r) 




I 


! 




• 


• 




• 


Y 




• LU • 


I 



We extend these cutting operations to forests u> <G OF by applying the B + operator to u> and then 
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cut it as a tree without using cuts of branches growing out of the root, before finally applying the 
B~ operator to R c (ui) to remove the added root. 



The coproduct of the Hopf algebra Hn (Section 2.3) can be formulated in terms of these 
cuts [28] . First one must extend the left admissible cuts to include the full cut of a tree, which 
cuts "below" the root, so that P^{t) is again r. The set of all left admissible cuts, including the 
full cut, is denoted by FLAC, and the coproduct A^r can be written as: 

A N (u) = P c {u)®R c {u). (29) 

ceFLAG 

Table [l] gives the result of this coproduct applied to all forests up to order 4. If we let An (to) 
consist only of forests resulting from not using the empty nor the full cut, we get A^(uj) = 
l(3uj + uj<g)l + Ajv(co). The operation Ajv is called the reduced coproduct. 

A formula for the substitution law. We will give a formula for the dual of the substitution, 
i.e. a formula for a+ , where (a*b,u>) — (b, a+ (uj)) . The formula is based on the pruning operation 
V on forests. 

Lemma 3.6 (Pruning). Letu andv be two forests. The dual of grafting, i.e. the operation defined 
by (v rv uj',oj) = (uj' ,V v (ui)), is given by: 

c€LAC{ui) 

The operation is called pruning. 

Proof. An elementary cut at an edge growing out of a node n of a forest uj is the dual operation 
of attaching trees via edges to the node n in a certain order, e.g. 



(oS) rx I, 



where the white nodes indicates where the attachment is done. The shuffling in P c (uj) corresponds 
to the dual of attaching forests in all possible ways to different nodes. Hence the dual of grafting 
is given by 

53 P c (u)<S>R c (uj). 

c<£LAC(lj) 



□ 



Theorem 3.7. We he 



**(*>)= E E ar(^(i))S + (anP c (^( 2 ))))a(iZ c (a; (2) )), 

(u)£i c c£LAC{lj (2) ) 



if uj ^ I. and (I) = I. Here A c denotes deconcatenation (Section 2.3). 

Note that using the magmatic product x defined in Section [2~3| this can also be written as: 

al = fj, o (fi x ® 7) o (oj ® a+ ® a) o (I ® A' N ) o A c , (30) 
where \x is concatenation, A^r is the coproduct in Hni and A' N (u) = Ajv(o;) — w<g)I for all forests 

UJ. 

Proof. We first prove the formula for ordered trees. Let 7Tot denote the projection of forests onto 
trees: 7tot(w) = X)t60t( t > uj ) uj - Recall that 

(o*a/,£j) = (uj' , (uj)) , (v rxu' ,u>) = (cj f , V v (u>)). 
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We have 

ir T{a+uj) = ^2 (t,o^uj)t = ^ (v rx •, aj(w)> v rx . 



tSOT ueOF 



(a-*- (z/ rx 1/ /^v • 
((a * ^) rx a, uj) v rx • 

i/eOF 

^ (a,V a * v w) v rxm 

i/£OF 

£ £ ( a * I /,P c H)( a ,fl c (a;)), 

c£LAC(uj) v£OF 



Hence, 



TroT(a^) = (",<i^M)(m.)a(ff( W )) 
= £ ((aHP c H)r>.)a(i? c (u;)) 

= £ i?+( a np c H)) a ( J R c H). 

ceL.4C(w) 

The general formula is established by the following calculation, where r is a tree: 

= X! ( v " r i a *( UJ )) ^ r 
= ^((a*z/)(a*r),u;)zAr 

(u)6A c i/,T 

51 (a*^(i))(7i"OT(a*W(2)))- 
(w)eA c 

As an example, the formula applied to the tree V yields 

e£(V) - a{V)B+(I) + a(*)a(l)B+(.) = a(Y). + a(.)a(I)I + a(«) 3 Y. 

See Table [2] where this formula is computed for all forests up to order 4, under the assumption 
that a is an infinitesimal character. 

Proposition 3.8. The map a* is a character for the shuffle product: 

a* (wi LU cu 2 ) = a*(wi) LU a* (w 2 ). 

Proof. The shuffle product is dual to the deshuffle coproduct, so the result follows from Lemma 
|3.4| by dualization. □ 



□ 



Remark 3.9. There is a similar formula for the substitution of B-series, only with the coproduct 
Ajv replaced by the Connes-Kreimer coproduct Ack = J2 c eAC(t) Pck( t ) ® ^cjf( r ) : 

a T Ar)= J2 B+(aT(pc K (T)))a(R c CK (r)) (31) 

ceAC(r) 



The proof of this formula is analogous to the proof of Theorem |3.7| This gives a recursive version 
of the coproduct in the substitution bialgebra Hcefm of [5] 
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3.3 Backward error analysis and modified vector fields 

Recall the results on backward error analysis in [7J: Given a B-series method Bf(a) there is a 
modified vector field f so that the B-series method applied to / generates the exact flow of /. 
Moreover, / can be written as a B-series with coefficients /3 satisfying /3 * 7Exact = ct, where 
7Exact is the coefficient function for the B-series of the exact flow, and * is the substitution law 
for characters in the Connes-Kreimer Hopf algebra Hqk- To generalize to LB-series, consider a 
numerical solution of the differential equation 

y' = f(y) ■ y (32) 

written in terms of a LB-series Bf(a) . We interpret it as the exact solution of a modified differential 
equation y' = f(y) ■ y. As in the classical case, it turns out that the modified vector field can be 
written as a LB-series / = Bf(f3). Furthermore, / is such that 

# e/ (M7Exact) = ( 33 ) 



where 7Exact represents the coefficients of the exact solution as described in Section |2.4| This 
result follows by applying Proposition |3.2| 



Theorem 3.10. Let Bf(a) be a LB-series method. There is a modified vector field f, given by 
f = £>/(/?) such that 

B}{lBxact)=B J {P). 

Moreover, 

(3 * 7 Exact = a. 

Example 3.11 (The exponential Euler method). The exponential Euler method is given by 

Vn+l = exp(hf(y n ))y n . 

In Example |2.17| the coefficients of the LB-series for this method was seen to be 7 = •. To get the 



backward error, we calculate j3 = Q(») o e, or log*(Q(»)) (cf. Section 2.4) 
In the classical setting this logarithm has been studied as log* (8), for a certain character S [6"Il3"0"]. 

4 Implementation 



As pointed out in Section |2.3| the set of forests T can be generated recursively using a magmatic 
product x defined on two forests u)\ and W2 by 

LJi X LU 2 = <JJ\B + {UJ2) (34) 

by starting with the empty tree I. Each forest in T can uniquely be written as a word in I and x . 
Recall that if u = u>i x uj 2 , then we call u)\ the left part, ujl, and LO2 the right part, ujr, of u). All 
the basic algebraic operations used to construct the substitution law can be formulated in terms 
of this product: 

Concatenation: wl = Iw = lo, and (u)\ x uj^) ^3 =ui x (W2W3). 

Shuffle: (jj[±\l = l\±\u) = ui, and uj\ LU 0J2 = LU L02L) x <^>ir + {^il LU ^2) x i^ir 

Coproduct: Ajv(I) = I® I, and A n (uj) = uj ® I + A n (uj l ) LU xAjv(wji) 



The formula ( 30 1 for the substitution law in Theorem |3.7| therefore lends itself well to implemen- 
tation. 
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Representing the free magma: One way to represent the free magma is by using well-formed 
words of parentheses '(' and ')'. A word w is well- formed if it is made of parentheses coming 
in pairs of one left and one right bracket, such that the left bracket appears on the left of the 
corresponding right bracket in w. For example, (())() is a well-formed word. The set of forests 
equipped with the product x is then isomorphic to this free magma via the recursion I = (), 

UJi X UJ 2 — (wi)W2- 

The authors have implemented a variant of the free magma, with elements represented by 
parentheses, and also the basic operations discussed in this paper. In future work, this implemen- 
tation will be used to do backward error analysis on interesting test cases, like the dynamics of 
rigid bodies. 
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A n (uj) 



1®I 

•(g)! + l<2>m 

I®I + + lol 
••(g)! + + 

i®i + «®i + 1®. + 1®! 

V®I + M<g>« + «®I + I® V 

•2(81 + 2—<gi» + •<E)l + «g>M 
L®I + t®» + + lot* 

•M(g)I + M(g>* + + I<8«« 

L)I + i(8« + I®I + »®1 + I® J 

+ + 2m<E>I + «®i + «<8 v + 1(8^ 

V(8l + Lg>» + + *2> V + I® V 
* S I / *«)I + *m(8« + + «(8 V + 1®*^ 

•2(81 + «I(8« + I»(8« + 2(8«« + 2«»®I + «(8«2 + »®I + 1(8 J 

L®I + 1(8» + 2(8«« + + I®L 

•V®I + 3«*«8« + M(g)M + 2«»®I + «(8«2 + »®V + I<8»V 
V»(8l + V(8» + + «<8)L + 1(8 V« 

n®i + l®« + «i(8« + 1®2 + 2M<g)M + «(8«2 + «®l + i®n 

•aLgl + 3m*(E)* + 2m®m + M(g)I + *®m« + «®«2 + I(E)m2 
•L®I + «I(8« + 2««(8«« + + + I(8«2« 

1*<E)I + + + *(g)«M + IiEiIm 

I + *M(8« + M0M + 



Table 1: Examples of the coproduct A 



N 
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i 

Y 
.1 

L 



Y 

J 
L 

.Y 
Y. 
II 
~I 
X 
!~ 



i 

«(•)• 

a(I). + a(.) 2 I 
a{*) 2 — 

a(i).+ 2a(.)a(I)I + a(.) 3 l 

a(Y). + a(«)a(!)l + a(.) 3 Y 

a(«I)« + a(»)a(I)«« + a(») 3 »I 

«(!•)• + a(«)a(I)«« + a(») 3 L 
a(») 3 «M 

a(l)« + 2a(»)a( )I + a(I) 2 I + 3a(») 2 a: 

ill 

a(Y). + a(«)a(I)I + a(«)a(Y)I + a(«) 2 a(I) ^ + ^ + a («) 4 Y 
a(^). + a(«)a(«!)I + a(«)a(!)I + a(«)a(Y)I + 3a(.) 2 a(I) Y + a(.) 4 ^ 

a(V). + a(.)a(I«)I + a(.)a(Y)I + a(I) 2 I + a(.) 2 a(I) ^ + ^ + a ('^ 
a(*V), + a(«)a(Y)I + a(«) 2 a(I) Y + a(.) 

a(i). + a(*)a(<!)! + a(«)a(!)M + 2a(«) 2 a(I)«I + a(«) 4 J 

a(I»)« + a(»)a(I«)I + a (•) «(!)•• + 2a(«) 2 a: (I)L + a(.)i 
a(.Y)« + a(.)a(«I)I + a(.)a(Y)«. + a(.) 2 a(I)«I + a(.) 4 .Y 
a(Y.). + a(«)a(L)I + a(«)a(Y)«. + a(.) 2 a(I)L + a(.) 4 Y. 
a(I) 2 .. + a(.) 2 a(I) (•! + !•) + a(.) 4 H 
«(«•!)• + a(«)a(«I)«« + a (•) 2 a (!)••• + a(») 4 »«I 
a(«2»)* + a (•) 2 a (!)••• + a(») 4 »I» 

a(L*)* + a(«)a(L)«« + a (•) 2 a (!)••• + a(») 4 l»« 

a(.) 4 - 



Table 2: Examples of the substitution character a+ , where a is an infinitesimal character, for all 
forests up to and including order four. 
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